function idat = fsb_remove_disc(idat,sandbox)

% FSB - DEV: Remove discontinuities from a 4D-Dataset
% This function belongs to fMRI Sandbox 
%
% EXAMPLE:
% idat = fsb_remove_disc(idat,sandbox)
%
% INPUT:
% idat:         4D image data set
% sandbox:      experiment information
%
% OUTPUT:
% idat:         a 4D dataset with the discontinuities removed
%
% CALLED BY:
% FSB.m
%
% NOTES:
% Does not yet work properly
% 
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate mean image and smooth data
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tmean = fsb_meancalc(idat,4);
midat = mean(idat,4);
s_idat = fsb_ssmooth_idat(midat,0.5);
s_idat(s_idat<1000) = 0;
s_idat(s_idat>1) = 1;
s_idat = single(s_idat);
idat = single(idat);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Now go over the data and look for discontinuities in the time course
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
h = waitbar(0,'Looking for discontinuities...');
for x = 1:max(sandbox.intrial(:,4))
    waitbar(x/max(sandbox.intrial(:,4)));
    trialseg = find(sandbox.intrial(:,4)==x);
    idatseg = idat(:,:,:,trialseg);
    idatseg = idatseg.*repmat(s_idat,[1 1 1 size(trialseg)]);
    idatstart{x} = mean(idatseg(:,:,:,1:3),4);
    idatstart{x}(idatstart{x}<500) = NaN;
    idatend{x} = mean(idatseg(:,:,:,end-3:end),4);
    idatend{x}(idatend{x}<500) = NaN;
    intrialdiff{x} = idatstart{x}-idatend{x};
    
    if x>1;
        intertrialdiff{x} = idatstart{x}-idatend{x-1};
    else
        intertrialdiff{x} = idatstart{x}-idatend{x};
    end

    idatseg = idat(:,:,:,trialseg)-(repmat(intertrialdiff{x},[1 1 1 size(trialseg)]));
    idatseg = idatseg.*repmat(s_idat,[1 1 1 size(trialseg)]);
    idat(:,:,:,trialseg) = idatseg;
    idatend{x} = mean(idatseg(:,:,:,end-3:end),4);
    idatend{x}(idatend{x}<500) = NaN;
    
end

close(h);
%idat(idat<500) = [];
idat = int16(idat);

tmean2 = fsb_meancalc(idat,4);
figure(701); cla;
set(gcf,'Name','Plot of mean volume intensity before (blue) and after (red) discontinuity removal');
plot(tmean); hold on; plot(tmean2,'r');

end
